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>>■ ABSTRACT 

^: 

We study evolution of various statistical quantities of smoothed cosmic density and 
velocity fields using N-body simulations. The parameter C = {V'^d'j /((V^) i^"^)) char- 
acterizes nonlinear coupling of these two fields and determines behavior of bulk velocity 
^ ■ dispersion as a function of local density contrast. It is found that this parameter de- 

. pends strongly on the smoothing scale even in quasi-linear regimes where the skewness 

in 
o 

o 
o 

c/3 . It is now commonly accepted that the large-scale structure in the universe is gravitation- 

^ ' ally evolved from small initial inhomogeneities. In addition to redshift surveys of galaxies and 

temperature anisotropics of the cosmic microwave background radiation, observational analysis of 
^ ' the cosmic velocity field would bring us fruitful information of our universe, such as, the density 

' parameter Jig or the matter power spectrum (Dekel 1994, Strauss & Willick 1995). In order to 

analyze the velocity field, however, we need to clarify a basic problem, whether the velocity field 
of galaxies is statistically different from that of dark-matter particles. The observed line-of-sight 
peculiar velocity field is not usually traced by the underlying density field but by galaxies, since the 
measurements of the distances are crucial element to determine the peculiar velocity field and are 
usually carried by using galaxies. This important problem in observational cosmology is broadly 
called the velocity bias. The velocity bias is often studied using semi-analytic models of galaxy 
formation (Cen &: Ostriker 1992, Kauffmann et al. 1999, Narayanan et al. 2000). But it is not 
easy to make definite predictions with this approach, because our understanding of astrophysical 
process related to galaxy formation reaches far from a quantitative level. 

The number density of galaxies is expected to be closely related to the density contrast of dark 
matter. Therefore we can obtain some insights about the velocity bias by analyzing relation between 
the velocity and density fields only of dark matter particles. This analysis is relatively simple, as 



parameter 5*3 = ((5^) / ^(5^) is nearly constant and close to the predicted value by the 
second-order perturbation theory. We also analyze weakly nonlinear effects caused by 
an adaptive smoothing known as the gather approach. 



1. INTRODUCTION 
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only gravitation is the relevant physical process for their evolution. When the initial fluctuations are 
isotropic random Gaussian distributed, these two fields at a same point are statistically independent 
in the framework of the linear perturbation theory. One of the author (Seto 2000a) studied weakly 
nonlinear evolution of the bulk (smoothed) velocity dispersion T,y(6) as a function of local smoothed 
density contrast S. He used the Eulerian second-order perturbation theory and the Edgeworth 
expansion method, and found that the constrained velocity dispersion 5]y((^) is written in the form 
T^yiS) (X {1 + C5 + 0{S'^)) with a parameter C ~ 0.2 ~ 0.3 for typical CDM models (see also 
Chodorowski &; Lokas 1997 for the velocity divergence field). 

However, numerical results given in Kepner et al. (1997) show that magnitude of bulk velocity 
is almost independent on the local density contrast S at nonlinear regimes (see their figs 2. a and 3. a). 
This behavior shows a remarkable difference from the above second-order prediction. Therefore we 
expect that the quantities C and Sy((^) show interesting behaviors by changing a spatial scale from 
linear to nonlinear, and would be useful measures to check performance of the perturbative analysis 
for weakly nonlinear evolution of the large-scale structure. In this article we study these quantities 
using N-body simulations and compare numerical results with analytical second-order predictions. 

This article is organized as follows. We begin by summarizing the second-order analysis of Seto 
(2000a) in §2.1. Smoothing operation is crucial for our numerical analysis and its proper treatment 
is very important to make a quantitative analysis. There are some variations for smoothing meth- 
ods. A mass-weighted smoothing scheme is convenient for analyzing N-body data and an adaptive 
smoothing scheme might be efficient to resolve cosmic fields, especially in sparse underdense re- 
gions. We perturbatively study these methods in §2.2 and §2.3. In §3 we derive an expression to 
estimate the sampling fluctuation expected in our numerical analysis. Numerical results are given 
in §4. Finally §5 is devoted to a summary. 



In this subsection we summarize the perturbative analysis of density and velocity fields and 
their correlations given in Seto (2000a). We denote the volume- weighted smoothing operator 
with radius R for a field f{x) as follows: 



where W{x, R) is a smoothing filter function. For simplicity we also use a notation fR{x) = [f{x)]ji. 
In this article we only employ the Gaussian filter defined as 



2. 



SECOND-ORDER PERTURBATION THEORY 



2.1. Smoothed Velocity Dispersion as a Function of Local Density 




(1) 




(2) 



and discuss velocity and density fields smoothed by this filter. 
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The smoothed (bulk) velocity dispersion for points x with a given smoothed overdensity 

^r{x) = 8 \s formally expressed as 

_ {VR{xf5n[5R{x)-5]) 

where is the Dirac's delta function and the bracket (■) represents an ensemble average. Let 

us evaluate this expression using the second-order Eulerian perturbation theory (Peebles 1980) and 
the Edgeworth expansion method (Matsubara 1994, Juszkiewicz et al. 1995, Bernardeau &: Kofman 
1995). We first expand the density and velocity fields as 

5{x) = 5i{x) + 52{x) + 6?,{x) + ■ ■ ■ , (4) 
V{x) = Vi{x) + V2{x) + V^{x) + ---, (5) 

where 5i{x) and Vi{x) arc the linear modes, 52{x) and V2{x) are the second-order modes, and 
so on. We can regard the order parameter of these expansions as the rms density fluctuation 
a = We assume that the primordial fluctuation is isotropic random-Gaussian distributed 

and the linear modes 5i [x] and Vi (x) obey this simple statistic. Then their probability distribution 
function P{Si,Vi) is determined by their covariance matrix {e.g. Bardeen et al. 1986). Although 
the linear velocity field Vi is given in terms of the linear density field 5i as Vi oc VA~^(5i, 
this matrix is orthogonal due to the statistical isotropy of fluctuations. Thus these four variables 
(density and three components of velocity fields) are statistically independent, as long as we discuss 
the density and velocity fields at a same point. Therefore at the linear order, velocity dispersion 
does not depend on the density contrast 6 and is given as 

^US) = {ViR-ViR) + 0{a^). (6) 

Next we evaluate the higher order correction using the multivariable Edgeworth expansion 
for the one point probability distribution function P{S, V) around its linear (gaussian) distribu- 
tion. After some tedious algebra, we obtain the following expression up to the first non-Gaussian 
correction 

^US) = {Vl)il + C5 + 0{a')), (7) 
where the coefficient C is defined by 

This coefficient C characterizes nonlinear couplings of the velocity and the density fields. Equation 
(7) shows that the first nonlinear correction for the velocity dispersion for points with a given 
density contrast S is simply proportional to CS. 

The leading order contributions for two factors (Vr ■ Vr) and {Sj^) in the denominator of C 
are given in terms of the matter power spectrum P{k) as follows 

m = j^^P{k)W{kRf + 0{a% (9) 
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= {Vr-Vr) = H^f I ^P{k)W{kRf + 0{a% (10) 



where H is the Hubble parameter and / is a function of cosmological parameters ilo and Aq, and 
well fitted by 

■,0.6 , -^0 A ^^0 



/-«r+^^i-^j, (11) 

in the ranges 0.03 < J^o < 2 and -5 < Aq < 5 (Lahav et al. 1991). The function W{kR) is 
a Fourier transformed filter with smoothing radius R. For the Gaussian filter (eq.[2]) we have 
W{kR) = exp(-/c2i?2/2). 

The ensemble average of an odd function of Gaussian variables (with vanishing means) leads 
to zero. We employ the second-order perturbation theory to evaluate the first nonvanishing con- 
tribution of the numerator (Vr ■ Vr5r). The leading order contribution for the numerator of C is 
written as 

{Vr ■ VrSr) = {ViR ■ VirS2r) + 2 {ViR ■ V2rSir) + 0(^6). (12) 

Using formulas for the second-order modes S2{x) and V2{x) {e.g. Fry 1984, Goroff et al. 1986), 
the above expression is written with the matter power spectrum P{k) as follows: 



u ( 5 u / k l\ 2 



fcZ 17 + 2 T + ^ +7" 



k + lu f3 u fk l\ An 



k{k^ + P + 2klu) [7 2 \ l kj 7 



+ Oia% (13) 



Here we have neglected extremely weak dependence on the cosmological parameters Qq and Aq that 
appears in the kernels of the second-order modes. We should notice that the parameter C = 0(1) 
does not depend on normalization of the power spectrum at its first nonvanishing contribution 
(see eqs.[9][10] and [13]). The factors [H fY cancel out between the velocity dispersion ay and the 
velocity-density moment {V r ■ V r6r). The parameter C is affected by the cosmological parameters 
mainly through its dependence on the power spectrum {e.g. the shape parameter F ~ fi/i for CDM 
transfer function). 

The parameter C characterizes the nonlinear mode couplings of velocity and density fields. 
It is interesting to compare this with a same kind of nonlinear quantity that is determined only 
by the density field 5. The skewness parameter ^3 = 0(1) is defined by the second- and third- 
order moments of density contrast 5 and written as (Peebles 1980, Pry 1984, Bouchet et al. 1992, 
Juszkiewicz, Bouchet & Colombi 1993, Bouchet et al. 1993, Bernardeau 1994, Baugh, Gaztanaga 
& Efstathou 1995, Kim & Strauss 1998) 



S3 = (14) 
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Here the variance (5^) is given in equation (9). For unsmoothed density field we have S3 = 34/7 
(Peebles 1980). After some algebra we can express the third-order moment ((J^) (Gauss filter) in 
terms of matter power spectrum P{k) as follows 

(5%) =6 f duj ^^^^Pik)Pil) exp[-(A;2 + p + klu)R^] 

The skewness parameter and its generalizations are very important to study various aspects of 
weakly nonlinear density field. With these parameters wc can discuss evolution of one-point PDF 
(Juszkiewicz et al. 1995, Bernardeau & Kofman 1995) or statistics of isodensity contour, such as, 
genus or area statistics (Matsubara 1994). It is known that the skewness parameter 5*3 obtained 
from numerical simulations shows a good agreement with second-order prediction up to regime 
(7^ ~ 1 {e.g. Hivon et al. 1995, Juszkiewicz et al. 1995, Baugh, Gaztanaga &; Efstathou 1995, 
Lokas et al. 1995). This fact is often used as a basis for reliability of second-order analysis. 

The skewness parameter has been also studied observationalily using various galaxy surveys 
(see Table 1 of Hui Sz Gaztaiiaga 1998). As we cannot observe the dark matter distribution S{x) 
directly, the relation between the number density fluctuation of galaxies Sg{x) and that of underling 
matter S{x) is crucially important to interpret the observed data. Their relation is called biasing 
and detailed theory of galaxy formation is required to understand it quantitatively. If the biasing 
relation is phenomenologically expressed as 

Sg{x) = h6{x) + ^j{S{xf - {S{xf)) + ■■■, (15) 

with numerical coefficients 6j, the first-order contribution of skewness parameter for the galaxy 
distribution Sg{x) becomes (Fry & Gaztaiiaga 1994) 

When we discuss velocity dispersion as a function of galaxy overdensity Sg , the coefficient C becomes 
(Seto 2000a) 

(''^'') = £ (17) 

Note that the above expression does not depend on the coefficient 62. 

In the following subsections we study various second-order effects of the moments S3 and C 
caused by smoothing operation. We denote results derived in this subsection by S^y and Cy with 
subscript V to specify the smoothing method which we employ, i.e., the simple volume- weighted 
smoothing method. 
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In Tables 1, 2 we present numerical values of Cy and Ssv for power-law models smoothed by 
Gaussian filter. The parameter S^y are given explicitly in terms of the Hyper geometric functions 
(Matsubara 1994, Lokas et al. 1995). For the top-hat filter we obtain simple expression as S^v = 
34/7 - (n + 3) (e.g. Bernardeau 1994). 



2.2. Second-Order Correction for the Mass- Weighted Smoothed Velocity Field 

In the previous subsection we have discussed the smoothed (bulk) velocity dispersion as a 
function of the smoothed density contrast S. Following standard procedure of the Eulerian pertur- 
bation theory, we have adopted the volume- weighted smoothing method (eq.[l]). However it is not 
a simple task to obtain the vohimc-wcightcd velocity field form N-body simulations due to their 
Lagrangian nature (Bernardeau & van de Weygart 1996, Kudlicki et al. 2000). In numerical anal- 
yses, a mass- weighted smoothing is a straightforward and convenient approach. In the followings 
we perturbatively investigate this approach (Seto 2000b for details). 

We consider particles' system mimicing N-body simulations. Particles are initially placed at 
grid points of the orthogonal (Lagrangian) coordinate system q. We assume that each particle has 
equal mass m{N) {N; number of particles in a simulation box). The simplest numerical method 
for calculating a smoothed velocity field (say at a point a;o) contains two steps. The first step is 
to sum up peculiar velocities of particles V{x{qj)) with weight W{x{qj) — Xo;R) determined by 
their distances to the point Xq in interest. The sum represents "momentum density" rather than 
"velocity". Thus the second step is to divide this sum with smoothed density at point Xq. We 
denote the final smoothed velocity field obtained in this manner by Vmass (xq ,R). It is expressed 
as follows: 

m{N)j:^W{x{g,)-xo;R) 

If we increase the total mimbcr of particles N with keeping their mean density p, we can replace the 
summation m{N) by three-dimensional integral with respect to the Lagrangian coordinate q. 
Next we change the coordinate system from Lagrangian q to Eulerian x and obtain the following 
result as 



lim m{N) V V{x{q,))W{x{q,) - Xo, R) = [ d\pW{x{q) - Xo; R)V{x{q)) (19) 

N—^00 ' J 

i 

= J d^xW{x-xo;R)p{l + d{x)}V{x) (20) 
= P{[V{xo)]r + [V{xo)S{xo)]r}, (21) 

where we have used the well-known fact that Jacobian \dq/dx\ of the coordinate transformation 
g — > a; is given by the local density 1 -|- S{x). As shown in the factor 1 -|- S{x) in equation (20), this 
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sampling is mass-weighted. In the same manner the denominator of equation (18) is given by 

N 

rn(N) 

N^oo 



lim m(Ar) V W{x{qi) - Xo, R) = p{l + [S{xo)]r). (22) 



In this case we can obtain the volume weighted quantity. With equations (20) and (22) the smoothed 
velocity field is expanded perturbatively as 

Vmassix, R) = [Vi{x)]r + [V2{x)]r + [V ,{x)5i{x)]r - [V r{x)]RMx)\R + 0{a\ (23) 

The linear term [Vi(a;)]R is unaffected by the smoothing method. The second term represents the 
volume- weighted smoothed velocity field at the second-order of perturbations. Let us extract out 
the new second-order correction terms V2m{x-,R) caused by the mass weighted smoothing as 

V2ra{x, R) = [Vi{x)5i{x)]r - [V i{x)]r[5i{x)\r. (24) 

Thus the velocity-density moment (V • V5) for smoothed fields Vmass ^'^d 5r is written as 

{Vlass^R) = {V15r) + 2 {ViR ■ V2m5m) . (25) 

The first term is given in the previous subsection up to the required order of a (see eq. [13]). The 
additional term caused by the mass- weighted smoothing is evaluated in terms of the power spectrum 
as 

2{[V]R.V2m6R) = 2i?V^/|^nW0{p + ^} 

X {exp (-(A;2 + l'^ + k ■ I)) - exp (-(fc^ + /2)) } , (26) 

For power-law models {P{k) oc A;") the correction term ACm caused by the above expression for 
the factor C is written with Hypergeometric functions (Seto 2000b) 

2{[V]R-V2n.SR) oz7/^^ + l " + 3 3 1\ il + n)^fn + 3 n + 3 5 1 



V2'2'2'4; 3 V2'2'2'4 



^ 2 ' 2 ' 2 , 
The values of this correction term for various n's are shown in table 1 



2.3. Adaptive Smoothing 

It has been pointed out that cosmic fields might be resolved more efficiently by using adaptive 
smoothing filters than traditional spatially fixed filters {e.g. Springel et al. 1998). Seto (2000b,c) 
studied second-order effects of this method for various moments of density and velocity fields. In the 
adaptive method the smoothing radius is determined by local density. We use a larger smoothing 
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radius at the sparse underdense region but use a smaller one at the high density region where 

nonlinear effects would be larger. Therefore, the perturbative treatment for the adaptive method 
might break down faster and its performance should be checked numerically. In this subsection we 
summarize second-order effects caused by the adaptive smoothing (Seto 2000b, c). Our analytical 
results are compared with numerical ones in §4. 

We first determine the adaptive smoothing radius so that number of particles within the radius 
becomes constant. This condition is written as 

R{xf[l + 6{x)]r(^x) =R^ = const, (28) 

where R is the standard smoothing scale. We can perturbatively solve this equation for the adaptive 
radius R{x) and obtain ^ 

R{x) = i? (l - H|)k + o{a')^ . (29) 

There are mainly two approaches for the adaptive smoothing, namely, the gather approach and 
the scatter approach (Hernquist & Katz 1989). Here we adopt the gather approach and basically 
calculate smoothed fields using operator ['J/j^aj)' After some algebra we obtain the second-order 
correction for the smoothed density field caused by the adaptive filter 

d2A{x,R) = -^5ir{x)^6Mx) + f ^ (Sfn) , (30) 

where the second term of the right-hand-side comes from the requirement {62A) = 0. This correction 
term gives new contribution AS^a to the skewness parameter 5*3 as 

3{52A{x,R)S,R{^r) 
- ^''^ 

= n + 3, (33) 

where we have used relation ((^x/j) 

R-in+S) for 

power-law models. As we can see from above 
derivation, the term AS^a does not depend on the filter functions. 

In the same manner the additional second-order terms for the smoothed velocity field is given 

as ^ 

V2a{x, R) = -^Sir{x)V^Vir{x). (34) 

We evaluate the new contribution ACa for the coefficient C caused by this correction term^ and 
obtain following result for power-law models as 



^In numerical analysis we use R{x) = R{1 + [5{x)]r) 
^Note that the contribution {S2AV1) vanishes. 
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3 dluR ^"^^^ 

(37) 



3 ' 

where we have used relation vahd for the Gaussian filter 

VV.«(., . (38) 

and a simple relation (T^ir) oc 

Since two second-order correction terms V2m and V2A are not coupled, their effects on mo- 
ments C are simply additive. For example, if we perform the mass-weighted smoothing using the 
adaptive filter, the moment C becomes 

C = Cv + ACm + AC A. (39) 
Same kind of additive relation holds for the skewness parameter. 



3. SAMPLE VARIANCE 

In the next section we numerically evaluate the factors C or S3 using N-body simulations. 
For each run of simulations we calculate the volume average A(Y,V) of a field Y(x) within the 
simulation box (cube) V as follows: 

A{Y,V)^^ J Y{x)d^x. (40) 

The ensemble average {A(Y,V)) of the estimator A(Y,V) coincides with the ensemble average 
{Y{x)) of the original field Y{x) and we have 

{A{Y,V)) = {Y{x)) . (41) 

But the measured value A{Y,V) fluctuates around its mean value {Y{x)) due to finiteness of the 
simulation box. This fluctuation can become important for analysis of numerical simulations. We 
call the root- man-square value D{Y,V) of this fluctuation as the sample variance 

DiY,Vf ^ {iAiY,V) - (Y))^) . (42) 

In this subsection we discuss the sample variances of various moments following Seto &; Yokoyama 
(2000). The sample variance is closely related to the number of statistically independent regions 
that is roughly determined by the correlation length of the field Y(x) and the size of the simulation 
box (or survey volume). Effects of the sample variance is generally more important for the velocity 
field than the density field, as the former is more weighted to large-scale fiuctuations and has a 
larger correlation length. 
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For fields Y{x) = V{x)'^d{x) and Y{x) = S{x)^ we have vanishing means (^d{x)V{x)^) = 
and (5(a;)^) = at the hnear order of perturbations as explained before^. However their sample 
variances do not vanish at the linear order and we have the following relation for relative fluctuations 
of their measured value A{Y,V) 

^^0(.-). (43) 

In contrast both the expectation values and the sample variances for the second moments Y[x) = 
5{x)'^ or Y{x) = V{x)'^ have contributions from the linear modes and we have 

^^0(1). (44, 

In the weakly nonlinear regime, the sample variance of the factor C = {V^6) /{{V^) would 
be dominated by fluctuations of the numerators (V^(5), as expected from equations (43) and (44). 

Here we evaluate the sample variance D{Y,V) caused by the linear modes for the moment 
Y{x) = V{x)'^6{x). We assume that our simulation box is a cube with side length L. After some 
algebra the sample variance can be written as 

2 77"4 

D{V'^6,Vf = -^^^P(fe)P(0P(|A; + Z|)iy(fei?)2M^(/i?)V(|fe + Z|)2 

x|^li-2y^|, (45) 



where the wave vector k is expressed as fe = 2Tr / L{nx,ny,nz) with integers Ux, riy and n^. Now 
let us approximate the above expression by replacing the summations with integrals, because it 
is not easy to evaluate these summations. We introduce the large-scale cut-off at wave-number 
kmin = 27r/L that is determined by the side length L of the simulation box V. Then we obtain 

2 TJ^ r r 

D{V^5,Vf ~ TTTTKT^ d^k dHP{k)P{l)P{\k + l\)W{kRfW{lRfW{\k + l\f 
(27r)"-L'^ Ju . J]^ . 



Due to the rotational symmetry around the origin this integral is simplified to a three-dimensional 
integral. In the same manner the sample variance for the third-order moments of density field is 
given as follows: 



^In this section we drop the subscript R (smoothing radius) and simply denote V{x) = [V{x)]r and S{x) = 
[5{x)U. 
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For the second moments (V(a;)^) or (^S{x)^) we have the following results 

Off 4: r If 

^ / d'kP{kfW{kRrk-\ D{5\ Vf ^ / d^kP{kfW{kRf- 

(48) 

In this section we made a crude estimation of the sample variance D{Y,V). We have only 
studied the fluctuations caused by the linear modes and used approximation of replacing summa- 
tions with integrals. Therefore our estimation is not quantitatively correct. However the present 
analysis would be useful to interpret numerical results given in the next section. 



4. NUMERICAL ANALYSIS 

We performed P'^M N-body simulations using 64^ particles in a cube with side length L = 1. 
The numerical code is provided by Couchman. We firstly investigated two scale-free models with 
spectral indexes n = 1 and in Einstein de-Sitter background and run several simulations for each 
models. We stopped calculations at the epoch when the matter fluctuation at the smoothing radius 
R = 1/80 {R: the smoothing radius for the Gaussian filter) becomes cr^ ci; 2 for the n = model 
and cj|j ~ 1 for the n = 1 model. Before describing numerical results, we mention two important 
effects that must be cared here. The first one is the sample variance described in the previous 
section and the second one is the artificial cut-offs induced in the initial fluctuations. 

We are interested in the velocity and density fields at weakly nonlinear regimes. As explained 
before, the sample variances of various statistical quantities become larger for a larger smoothing 
radius. In figure 1 we estimated the relative fluctuations D{Y,V)/ {Y) caused by the linear modes 
using expressions given in §3 at the final epoch of each simulation. Apparently, the moment 
^V^(5) has the largest fluctuations. The fluctuation can be as much as the expectation value itself, 
i.e., 100% variance for smaller a^{R). The sample variance is very sensitive to the large-scale 
fluctuations. Roughly speaking, the sample variance is more important for n = model than for 
n = 1 model, as the former has more large-scale powers. Because of the same reason, the velocity 
dispersion {V'^) has larger fluctuations than the variance of density fluctuations {S'^)- 

There are two cut-off scales for initial matter fluctuations in our simulations. One is the 
Nyquist frequency kNyq = ttN/L determined by the initial separation of particles and the other one 
is the wave number kmin determined by the side length of the simulation box as kmin = 27r/L. Our 
analytical predictions given in §2 might be different from numerical results due to these artificial 
cut-offs. To estimate their efl'ccts we calculate the factors Cy and ACm with a modified power 
spectra Peff{k) mimicing N-body simulations as follows {e.g. Seto 1999): 

{0 {0<k< kmin) 
A;" {km,n <k< kNyq) (49) 
{kNyq < k). 
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Analysis for the n = — 1 model is particularly interesting as we have Cy = ACm = AC^ = and 
velocity dispersion does not depend on the density contrast even at the second-order of perturbation. 
However, the wave- number integral (eq.[10]) for the velocity dispersion diverges in the A; ^ limit 
and the velocity field would be strongly affected by the artificial large-scale cut-off at kmin- Actually 
we have confirmed numerically that the factors Cy and ACm for this model with the modified 
spectrum Peff{k) would considerably deviate from the original scale-free model. We apply this 
analysis for models with indexes n = and 1, and found that the factors Cy and ACm could 
deviate largely from predictions for the original scale- free models at R < 1/80 or > 1/20. 
Thus we limit our analysis only for smoothing radius with 1/80 < R < 1/20. 

We have analyzed N-body data at the final epoch for each run. As our simulations are basically 
self similar, we can expect that physics is characterized only by the magnitude of nonlinearity 
CT^{R) = {^r}- Thus we use the variance <t^{R) rather than the radius R to represent the scale 
dependence. For n = model we have performed three runs. Our numerical results are given in 
figure 2. In the upper left panel the factors Ssv and C are shown as a function of nonlinearity 
a'^{R). The symbols represent the average values of three runs and the error bars are the variances 
of their fluctuations. The magnitude of error bars is roughly understood by checking the sample 
variances caused by the linear modes given in figure 1 where we have added numerical results from 
simulations for the quantity (V'^S'j. 

As shown in literatures (e.g. Hivon et al. 1995, Juszkiewicz et al. 1995, Lokas et al. 1995), 
the skewness parameter S^v at ct^{R) < 1 is nearly constant for n = 1 model and close to the 
predictions of the second-order perturbation theory S^y = 3.14. In contrast the factor C depends 
strongly on nonlinearity cr^(i?) and we have C ~ 0.05 at (t^{R) — 1. This value is only ~ 40% of 
the second-order prediction C = 0.12. The numerical results for C become close to the analytical 
prediction at cr'^{R) = 0.2 where the sample variance is fairly large. In figure 2 we also show the 
constrained velocity dispersion T,'y{5)/ay in the range of density contrast —1.5a{R) < 5 < 1.5a(R). 
We have made a same kind of averaging operation for three runs as in the case of the factors S^y and 
C. The dashed lines are the analytical predictions given in equation (7) whose slope is C 2± 0.12. 
The numerical results is close to the analytical prediction at small cr^{R), but the slope of numerical 
results become smaller for larger cr{R) as expected form behavior of the factor C. 

Let us show same analyses for the model with n = 1. For this model we performed seven 
runs. Results are shown in figure 3. The measured parameter Ssv is nearly constant but somewhat 

(~ 10%) smaller than the analytical prediction S^v = 3.03. The factor C is again a decreasing 
function of nonlinearity a'^{R). However its convergence to the analytical prediction C = 0.2 seems 
to be slower. Even at cr^(-R) = 0.04 where the sample variance is considerably large, the measured 
value is much smaller than the analytical one. The constrained velocity dispersion Sy((^) behaves 
in the same manner as the model with n = 0. The effective slope becomes smaller for larger 
cr^(i?). The difference between the second-order prediction and the numerical results is caused 
by the nonlinear effects that cannot be traced by the second-order perturbation. For given o"^(i?) 
this difference is larger for the model with n = 1 than the one with n = 0. This is reasonable as 
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the former has more small scale power initially and nonlinear effects would be larger. Nagamine, 
Ostrikcr &: Cen (2000) have numerically studied the cosmic Mach number (Ostriker & Suto 1990) 
as a function of overdensity using an LCDM hydrodynamical simulation. They have provided two- 
dimensional contour plots for the probability distribution functions of the smoothed dark-matter 
overdensity and the bulk velocity \Vii\ at various smoothing scales R = 5,10 and 20/i~^Mpc 
(see their figures 10-12). Comparing these figures we can observe the scale dependence of Sy(5) 
discussed so far. 

We have also performed N-body simulations (5 runs) for a CDM model. Our model parameters 
are fixed at JIq = 1-0, h = 0.5 ^ and the CDM shape parameter F = 0.5. We used the CDM 
transfer function of Bardeen et al. (1986) with the primordially n = 1 spectrum and normalized 
the linear spectrum by as = 0.8 (erg is the rms density fluctuation within a top-hat sphere of 
radius 8/i~^Mpc). As the velocity field is heavily weighted to large-scale fluctuations, we have 
adopted a L = 400/t~^Mpc simulation box. The prediction C by the second-order theory at the 
weakly nonlinear regime is C ~ 0.13 and close to the previous result for the n = model. Our 
numerical results arc C = 0.094 ±0.06 for R = 10/i~^Mpc {a'^{R) = 0.07), and C = 0.068 ±0.01 for 
R = 5/i~^Mpc {a'^{R) = 0.37), which are consistent with the prediction of the second-order theory. 
The mean values and the error bars depend on nonlinearity as in the case of the n = model. 

Hui & Gaztahaga (1999) have pointed out that a nonlinear combination / of estimators A{Yi, V) 
(see eq.[40]) should be analyzed with sufficient care. We have following inequality even in the case 
of the unbiased estimators {A{Yi,V)) = (Yi) : 

/((Yi) , ■ ■ ■ , (y„)) = /((^(Yi, V)) , • • • , (A(F„, V))) ^ {f{A{Y^, V), • • • , A(y„, V))) . (50) 

This inequality leads to a biased estimation for the nonlinear function /((l^)). When the variance 
^ of density contrast smoothed in survey volume V is much smaller than unity, this bias ASsv for 
the estimation of the skewness Ssv is approximately given as follows (Hui & Gaztahaga 1999) 

where the coefficients {ui} are determined by the power spectrum. As we use data in full simulation 
boxes, there are no mean density fluctuations ^ = and the above expansion (51) vanishes. There- 
fore we evaluate the estimation bias by directly comparing numerical values f{{A{Yi,V)) , ■ • • , (^(1^, V))) 
and {f{A{Yi,V), ■ ■ ■ ,A{Yn,V))) (the ensemble average (■) is taken for different realizations of sim- 
ulations). Wc find that differences between these two arc very small both for Sy and Cy- For data 
presented in figures 2 and 3 the difference is smaller than 1%. This is due to the fact that the 
fluctuations of their denominators (expressed by ((^^) and {V'^)) are very small as expected from 
figure 1. This estimation bias would be more important in actual observational situations. 



''As commented in §2.1, our results for the parameter C would depends weakly on the background cosmological 
parameters. 



(51) 
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In figure 4 we present numerical results for the factors S3 and C that were evaluated using 
the adaptive filter described in §2.3. As shown in the left panels, the skewness parameter 
shows reasonable agreement with the second-order perturbation theory. But the parameter C 
again shows strong dependence on nonlinearity cr^{R) even in this regime and its convergence to 
the analytical prediction is slow. These behaviors of S and C are similar to the previous case 
of non-adaptive smoothing. The parameter C for the model with n = becomes close to the 
second-order result C = 0.45. For n = 1 model, the second-order correction ACa = 0.67 caused by 
the adaptive smoothing seems not to match the numerical result even at the scale ~ 0.04. In 
the adaptive method we use a smaller smoothing radius at high density regions. Thus numerical 
results are expected to be more contaminated by strongly nonlinear effects than the previous simple 
smoothing method especially for n = 1 model. The numerical result is estimated as ACa — 0.44 
from Cv + ACm + ACa = 0.56 (figure 4) and Cy + ACm = 0.12 (figure 3). 

Finally we comment on possibility of the determination of the parameter Cy by the actual 
observations. As only the line-of-sight velocity field V\\ {x) is observable, we evaluate the fluctuations 

of linear modes for the moment (v^^5^. For this evaluation we cannot simply apply results derived 
under the periodic boundary condition in §3. We estimate the sample variance D(V|^(5, V) using a 
method similar to Seto & Yokoyama (2000). As an example we have examined a flat cold-dark- 
matter model with cosmological parameters h = 0.7, f^o = 0.3 and Aq = 0.7. Normalization of 
the power spectrum is determined by abundance of rich clusters as erg = 0.501]o°-^^+°-^2"" (Eke, 
Cole &: Frenk 1996). The Gaussian smoothing radius is fixed at i? = 12/i^^Mpc. The second-order 
prediction of the parameter C for the primordially scale-invariant CDM power spectrum in a weakly 
nonlinear regime is close to that for the n = model. Therefore we assume 

V?d) = ^— ^ ~ J - 0-04 (V^) (6^) . (52) 



Using these relations we obtain D{V^^S,V)/ ^^^<^/ ~ 8 for a spherical volume with survey depth 
80/i~^Mpc. Thus we can hardly measure the parameter C from current catalogs of peculiar ve- 
locity field, such as, Mark III (Willick et al. 1997) or SFI (Haynes et al. 1998). The relative 
fluctuation becomes close to unity only at a survey volume with depth ~ 300/i~^Mpc. In the ac- 
tual observational situations the effects caused by the nonlinear estimators (eq.[50]) would be also 
important. 



5. SUMMARY 

We have investigated various statistical quantities of smoothed velocity and density fields. 
Analytical predictions based on second-order perturbation theory are compared with numerical 
results obtained from N-body simulations. Our primary target is the velocity dispersion T,y{6) as 
a function of local density contrast S at weakly nonlinear scales. At the linear order with isotropic 
random Gaussian initial condition the dispersion Sy(5) does not depend on S. However second- 
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order modes generate a correction term proportional to S with its slope C = (V^d) /((V^) (^^)). 
We have confirmed this dependence at small cr{R) and found that the parameter C depends largely 
on nonlinearity (j{R) even in the regime where the skewness parameter is nearly constant value 
and close to prediction by second-order perturbation theory. 

We have also studied second-order effects caused by various smoothing methods. Proper 
treatment of smoothing is crucial to quantitatively discuss weakly nonlinear effects of cosmic fields. 
A mass weighted smoothing is convenient to analyze velocity field traced by particles as in N-body 
simulations. We have basically used this method for our numerical analysis. We found that the 
parameter C actually approaches to the value that includes effects of mass-weighted smoothing 
in the limit of small a(R). But a large smoothing radius is required to recover the second-order 
result. For example we have to take the radius R with < 0.1 for the n = model and 

CT^{R) < 0.01 for the n = 1 model. Adaptive smoothing methods might be useful to resolve 
cosmic field especially in underdense regions. The skewness parameter ^3 obtained numerically 
with an adaptive method (gather approach) is closed to the corresponding second-order prediction 
up to cr^{R) — 1. The skewness parameter and its generalizations are closely related to weakly 
nonlinear evolution of genus or area statistics of isodensity contour (Matsubara 1994) for which the 
adaptive method would be effective (Springel et al. 1998). Our result is encouraging for application 
of second-order analysis for density field smoothed by the adaptive filter. 

It has been known that the second-order predictions at (7{R) ~ 1 are more accurate for the 
skewness of the real space density field than that of the velocity divergence field (xV-V ( Juszkiewicz 
et al. 1995) or the redshift space density field (Hivon et al. 1995) . But we should notice that calcula- 
tion of the coefficient C = (V^S) /((V^) {S'^)) directly requires information of the peculiar velocity 
field as a vector field. We found that performance of the second-order prediction is generally worse 
for the factor C than for the skewness 6*3. This fact seems interesting. As the density field is more 
weighted to the small-scale fluctuations than the velocity field, we can naively expect that highly 
nonlinear effects would be stronger for the density field. Therefore our numerical investigation 
provides applications of the second-order perturbation theory with an important caveat. 

We thank Naoki Yoshida for discussion about numerical schemes and an anonymous referee for 
valuable comments to improve this manuscript. This work is supported by Japanese Grant-in- Aid 
for Science Research Fund of the Ministry of Education, Science, Sports and Culture Grant Nos. 
0001461 and 11640235. 
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TABLE 1 

VELOCITY-DENSITY MOMENT C (GAUSSIAN FILTER) 



spectral index n 


-1 





1 


2 


Cv 





0.24 


0.31 


0.38 


ACm 





-0.12 


-0.11 


0.053 


ACa 





0.33 


0.67 


1.0 


Cv + ACm 





0.12 


0.20 


0.43 


Cv + ACm + ACa 





0.45 


0.87 


1.43 




TABLE 2 






SKEWNESS Ss (Gaussian filter) 




spectral index n 


-1 





1 


2 


Ssv 


3.5 


3.1 


3.0 


3.1 


AS3A 


2.0 


3.0 


4.0 


5.0 


Ssv + AS3A 


5.5 


6.1 


7.01 


8.1 
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CJ2(R) 



Fig. 1. - Relative fluctuations caused by the linear modes in our simulations (see §3 for their 
estimation). The solid lines correspond to the field Y{x) = V{x)'^5{x), the dotted to Y{x) = d{x)^, 
the short-dashed to Y{x) = V{x)'^ and the long-dashed to Y{x) = S{x)'^. We numerically estimate 
the fluctuations D{V'^S,V) from different realizations of simulation (see figures 2 and 3) and plot 
its relative magnitude with respect to the expectation value (V'^S) (predicted by the second-order 
perturbation theory) by the filled squares. 
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Fig. 2. — Two factors C, S^v and the constrained bulk velocity dispersion S^(^) for n = model. 
(1) Upper-left panel: The mean values (filled symbols) and error bars are obtained from three runs 
of simulations. The squares represent the factor C in the form of 20C, and the triangles the skewness 
S3V Predications by the second-order perturbation theory are shown by the dashed lines (C = 0.12 
and S3V = 3.14). (2) Other three panels: The constrained velocity dispersion is presented using 
the ratio 'Ey{S)/ay in the range —1.5a{R) < S < 1.5a{R). The short-dashed lines are predictions 
by the Edgeworth expansion method with the factor C from second-order perturbation theory, and 
the long-dashed lines are drawn with C obtained from numerical simulations. 
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Fig. 3. — Numerical results for n = 1 model. Figures are given in the same manner as figure 2. 
Second-order perturbation predicts C = 0.20 and Szv = 3.03. 
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Fig. 4. — Effects of the adaptive smoothing for moments ^3 and C. We present the averaged values 
as in Figures 2 and 3. The upper panels correspond to n = model and lower panels to n = 1 
model. The dashed lines represent the second-order predictions: Sz = 6.1, C = 0.45 for n = 
model and S3 = 7.0, C = 0.87 forn = 1 model. 



